### This script will receive an argument from a .sh file, which is year
### All BRT ASCII's for that year will be read in and stored in a single dataframe
### That data frame will be saved as a .csv file
### This is an intermediate step in the extreme weather analysis

library('SDMTools')

#### Location of BRT tmax ASCII files 

tmax.dir = '/home/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/brtpredsfinal/maxgzip/'

### Need to create base.pos with all positions to intersect with their respective lat/longs
### First read in the DEM ASCII as the template

base.asc = read.asc('/home/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/dem_WTplusbuffer_LatLong_WGS1984_250mres.asc')

### Identify all positions where data exists (n should equal 499085)

base.pos = as.data.frame(which(is.finite(base.asc), arr.ind=T))

### Convert all remaining finite data in base.asc to NA

base.asc[,] = NA

### Get lat longs for data extraction

base.pos$lat = getXYcoords(base.asc)$x[base.pos$row]
base.pos$long = getXYcoords(base.asc)$y[base.pos$col]

### List files to intersect

tfiles = list.files(paste(tmax.dir,sep=''),pattern='.asc.gz',full.names=T,recursive=T)

### Create an object 'outdata' to bind data onto

outdata = base.pos

### Position tracker

i=1

for (yearx in c(1951:2011))

	{
	
	### Subset files 
	
	ttfiles = tfiles[which(substr(tfiles,nchar(tfiles)-14,nchar(tfiles)-11)==yearx)]
	
	for (tfile in ttfiles)
	
		{
		
		### Read in an ASCII
		
		tasc = read.asc.gz(tfile)
		
		### Round of ASCII data to two decimal points
		
		tasc = round(tasc,2)
		
		### Extract data from points in base.pos
		
		tdata = extract.data(cbind(base.pos[,3],base.pos[,4]),tasc)
		
		### Bind to outdata
				
		outdata = cbind(outdata,tdata)
		
		#### Change the name of the newly bound column

		names(outdata)[i+4]=paste(substr(tfile,nchar(tfile)-14,nchar(tfile)-7),sep='')
		
		### Change position tracker
		
		i=i+1
		
		#### Report progress
		
		cat('\n',i,'\n')
		
		}
		
	### Write outdata to the appropriate directory

	write.csv(outdata,file=paste(tmax.dir,yearx,'/allpos_alldays_maxtemp_',yearx,'.csv',sep=''), row.names=F)	
	
	### Clear outdata
	
	outdata = base.pos
	
	### Reset position tracker
	
	i=1
		
	}
	
### Close loop
		
			
			
			
			
			
			
			
			
			
			
			
			
			
			